!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Common layout of the time integrators
!!
!! \n
!!
!! Loosely speaking, the generic ODE problem is described by the
!! abstract type <tt>c_ode</tt>, while the generic ODE integrator is
!! described by the abstract type <tt>c_tint</tt>. To define a
!! specific ODE, the user extends \c c_ode, while to provide a
!! specific time integrator a library extends \c c_tint. Hence, \c
!! c_ode represents a continuous problem, while \c c_tint represents a
!! discretization method.
!!
!! \section sectODE The continuous problem
!!
!! Representing the generic ODE as
!! \f{displaymath}{
!!  y' = f(t,y),
!! \f}
!! it is clear that in principle we need two types: \c c_stv
!! corresponds to \f$y\f$ while \c c_ode corresponds to the function
!! \f$f\f$. Often, however, some scratch area is required to compute
!! \f$f(t,y)\f$, which complicates the overall organization. This
!! scratch area should not be included in \c c_stv, since it does not
!! need to be replicated for different time levels, and it should also
!! not be included in \c c_ode because \c c_ode objects are not
!! modified during the time integration (more precisely: they are not
!! modified in the computation of the rhs, while the work arrays can
!! be modified during the computation of the time step). Therefore, a
!! third type is required: \c c_ods (ODE scratch). While a detailed
!! interface is specified for \c c_stv and \c c_ode, the abstract \c
!! c_ods type is empty, since it completely depends on the specific
!! equation. The only thing that is ensured by the general layout is
!! that such an object is passed to the functions contained in \c
!! c_ode.
!!
!! In a typical problem, one would distribute the problem data as
!! follows:
!! <ul>
!!  <li> <tt>c_stv</tt>: each variable having a time evolution, in
!!  particular all the prognostic variables
!!  <li> <tt>c_ode</tt>: everything defining the ODE and remaining
!!  constant during the integration, such as  the right-hand-side of
!!  the ODE (in the dedicated fields of the abstract \c c_ode class),
!!  as well as any parameter characterizing the ODE
!!  <li> <tt>c_ods</tt>: anything that can be assimilated to a
!!  diagnostic variable, which could, at least in principle, be
!!  eliminated from the problem.
!!  <li> when a PDE is discretized with the method of lines, there are
!!  two equally sensible choices for the finite element grid and
!!  basis: they can be included either in \c c_stv or in \c c_ode;
!!  <ul>
!!   <li> when they are included in \c c_stv, they should be included
!!   as pointers, so that different \c c_stv objects share the same
!!   grid and basis without duplications
!!   <li> when they are included in \c c_ode, they can be included
!!   either as pointers or as components, since in general there is
!!   only one \c c_ode object
!!   <li> a possible solution is including them as pointers in both \c
!!   c_stv and \c c_ode: then, for instance, <tt>c_stv\%scal</tt> will
!!   use <tt>c_stv\%grid</tt> and <tt>c_ode\%rhs</tt> will use
!!   <tt>c_ode\%grid</tt>, both of them actually using the same grid.
!!  </ul>
!! </ul>
!! \warning During a call to <tt>c_tint\%step</tt> both \c c_stv and
!! \c c_ods are modified, however only the values contained in \c
!! c_stv are meaningful and can be associated with a specific time
!! level, while the values of \c c_ods should be ignored. If they are
!! required, for instance because some output for the corresponding
!! diagnostic variables is required, then they must be computed
!! explicitly by the user using \c c_stv; this operation anyway does
!! not affect the time integration.
!!
!! The right-hand-side \f$f\f$ is a function, and thus should require
!! one corresponding method. In fact, \c c_ode includes two methods,
!! corresponding to the "viewpoints" of explicit and implicit time
!! integrators. Often, each solver uses only one of these two methods.
!!
!! \section sectTI The numerical method
!!
!! The fields included in \c c_tint should suffice for a very large
!! class of time integrators; notice that a particular time integrator
!! does not necessarily use all them.
!!
!! Each time integrator must be initialized once, then it can be
!! called an arbitrary number of times to advance the solution, and
!! finally it must be cleaned up.
!!
!! The initialization subroutine can be very simple or rather complex,
!! depending on the method. Typically, for one-step methods the
!! initialization boils down to allocating the work arrays. For
!! multi-step methods however the initialization includes advancing
!! the numerical solution for the first time steps: this motivates the
!! rather rich interface of \c i_ti_init (see also the fields
!! <tt>t_ti_step_diag\%bootstrap_steps</tt>).
!!
!! \section details Implementation details
!!
!! Many interfaces specify <tt>intent(inout)</tt> instead of the more
!! natural <tt>intent(out)</tt> for the \c c_stv, \c c_ode and \c
!! c_ods arguments because the specifications of such objects will
!! likely include large allocatable components, and it would be
!! inefficient to reallocate them at each time step.
!!
!! For the same reason, we also avoid functions in favour of
!! subroutines, because functions would require frequent
!! allocations/deallocations of the function results, while the
!! subroutine arguments can be allocated once and then passed around
!! thanks to the \c inout intent.
!!
!! The two types \c t_ti_init_data and \c t_ti_step_diag provide a
!! general interface to pass parameters to and from the time
!! integrators. The aim is simplifying switching among different time
!! integrators, although some specific, integrator depending settings
!! can not be eliminated.
!<----------------------------------------------------------------------
module mod_time_integrators_base

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   mpi_comm_world

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_time_integrators_base_constructor, &
   mod_time_integrators_base_destructor,  &
   mod_time_integrators_base_initialized, &
   c_ode, c_ods, c_tint, &
   t_ti_init_data, t_ti_step_diag

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> ODE problem
 !!
 !! \note Some procedures (namely \c solve) are required only by some
 !! time integrators. We provide here a dummy version of such
 !! subroutines, so one doesn't have to define them as long as he/she
 !! doesn't use the time integrators requiring them.
 type, abstract :: c_ode
 contains
  !> Evaluate the right hand side \f$f(t,y)\f$
  !!
  !! The optional argument \c term is introduced to allow the
  !! implementation of splitting methods. For these methods, the right
  !! hand side is decomposed as
  !! \f{displaymath}{
  !!  f(t,y) = \sum_{i=1}^s f_i(t,y)
  !! \f}
  !! and the argument \c term is used to select
  !! <ul>
  !!  <li> <tt>term(1) .eq. 0</tt> \f$\mapsto\f$ \f$f(t,y)\f$
  !!  <li> <tt>term(1) .gt. 0</tt> \f$\mapsto\f$ \f$f_{{\tt
  !!    term(1)}}(t,y)\f$.
  !! </ul>
  !! When this argument is absent, the whole rhs is implied. The
  !! second element <tt>term(2)</tt> is always equal to \f$s\f$, the
  !! number of terms required by the time integrator. This is useful
  !! to implement \c rhs in such a way that it can work both with a
  !! three-term splitting and, for instance, a two-term one, where
  !! some terms are grouped together.
  procedure(i_rhs), deferred, pass(ode) :: rhs
  !> Solution of the implicit problem \f$ x - \sigma f(t,x) = b \f$.
  !!
  !! This subroutine is required only for implicit solvers, in which
  !! case it solves the implicit problem for arbitrary
  !! \f$\sigma\in\mathbb{R}\f$ and right hand side \f$b\f$. An
  !! additional argument \c xl is provided to pass a state which can
  !! be used to approximate \f$f\f$ (typically by linearization), thus
  !! solving
  !! \f{displaymath}{
  !!  x - \sigma \tilde{f}_{x_l}(t,x) = b
  !! \f}
  !! where \f$\tilde{f}_{x_l}\f$ approximates \f$f\f$ in a
  !! neighbourhood of \f$x_l\f$.
  !!
  !! \note Sometimes the solution of the implicit problem requires an
  !! inital guess. In such a case, all the time integrators provide
  !! two possible values:
  !! <ul>
  !!  <li> the extrapolation \f$x_l\f$, computed by the method itself,
  !!  <li> the value of \c unew passed in the call to the
  !!  <tt>interface(i_ti_step)</tt> subroutine (or the value \c u0
  !!  passed in the call to the <tt>interface(i_ti_init)</tt>
  !!  subroutine for the initialization phase).
  !! </ul>
  !! The first option is useful when \f$x_l\f$ is a good approximation
  !! of \f$x\f$, and relies on the implementation of the
  !! time-integrator (which typically uses high-order extrapolation).
  !! The second option gives the user the highest control on the
  !! definition of the initial guess, for example by using the
  !! solution at the previous time level.
  procedure, pass(ode) :: solve
 end type c_ode

 !> ODE scratch data (diagnostic variables)
 !!
 !! This type is not abstract because for simple problems an empty
 !! type could suffice, hence this type can be used directly. Clearly,
 !! whenever some scratch memory <em>is</em> required, this type
 !! should be extended.
 type :: c_ods
  ! add your fields as required
 end type c_ods

 !> ODE integrator
 type, abstract :: c_tint
  real(wp) :: dt !< time step
  !> Storage available to the time integrators to deal with
  !! temporaries and/or multiple time levels.
  !!
  !! This array will be allocated to a size which depends on the
  !! requirements of the specific time integrator and to a type which
  !! is taken from the problem unknown.
  class(c_stv), allocatable :: work(:)
 contains
  procedure(i_init),  deferred, pass(ti) :: init
  procedure(i_step),  deferred, pass(ti) :: step
  procedure(i_clean), deferred, pass(ti) :: clean
 end type c_tint

 !> Some useful parameters that <em>could</em> be used by various time
 !! integrators. Notice that the field names are rather generic, each
 !! time integrator will use these values in a specific way (or will
 !! ignore them).
 type t_ti_init_data
  integer  :: dim = 10
  real(wp) :: tol = 1.0e-6_wp
  !> A collection of integer input parameters.
  !!
  !! They have a different meaning for different time integrators.
  integer, allocatable :: ii1(:)
  !> A collection of real input parameters.
  !!
  !! They have a different meaning for different time integrators.
  real(wp), allocatable :: ir1(:)
  integer  :: mpi_id = 0
  integer  :: mpi_comm = mpi_comm_world
 end type t_ti_init_data

 !> This type can be used to retrieve various diagnostics from the
 !! time integrator step.
 !!
 !! Some general purpose fields are introduced here; one can consider
 !! extending this type to provide additional diagnostics for specific
 !! time integrators. To deal with allocatable components, objects of
 !! this type should be passed to the time integrators initializer
 !! with intent out and then to the time integrator step with intent
 !! inout.
 !!
 !! Each time integrator can set the fields <tt>name_d*</tt> according
 !! to the returned diagnostics.
 type t_ti_step_diag
  !> This field is used in the initialization of multi-step methods to
  !! indicate how many steps are performed in the start-up phase.
  integer :: bootstrap_steps = 0
  logical :: max_iter
  integer :: iterations
  character(len=100) :: name_d1 = 'd1'
  character(len=100) :: name_d2 = 'd2'
  character(len=100) :: name_d3 = 'd3'
  real(wp), allocatable :: d1(:)
  real(wp), allocatable :: d2(:,:)
  real(wp), allocatable :: d3(:,:,:)
 end type t_ti_step_diag

 !> tendency computation
 abstract interface
  subroutine i_rhs(tnd,ode,t,uuu,ods,term)
   import :: wp, c_ode, c_stv, c_ods
   implicit none
   class(c_ode), intent(in)    :: ode !< ODE problem
   real(wp),     intent(in)    :: t   !< time level
   class(c_stv), intent(in)    :: uuu !< present state
   class(c_ods), intent(inout) :: ods !< scratch (diagnostic vars.)
   class(c_stv), intent(inout) :: tnd !< tendency
   !> term to be evaluated and number of terms in the integrator
   integer,      intent(in), optional :: term(2)
  end subroutine i_rhs
 end interface

 !> Generic init subroutine of the time integrators
 abstract interface
  subroutine i_init(ti,ode , dt,t,u0,ods , init_data,step_diag)
   import :: wp, c_tint, c_ode, c_stv, c_ods, &
             t_ti_init_data, t_ti_step_diag
   implicit none
   class(c_ode),  intent(in)    :: ode  !< ODE problem
   real(wp),      intent(in)    :: dt, t !< time step, initial time
   class(c_stv),  intent(inout) :: u0   !< initial condition
   class(c_ods),  intent(inout) :: ods  !< scratch (diagnostic vars.)
   class(c_tint), intent(out)   :: ti   !< time integrator
   type(t_ti_init_data), optional, intent(in)  :: init_data
   type(t_ti_step_diag), optional, intent(out) :: step_diag
  end subroutine i_init
 end interface
 !> Generic step subroutine of the time integrators
 abstract interface
  subroutine i_step(unew , ti,ode,t,unow,ods , step_diag)
   import :: wp, c_tint, c_ode, c_stv, c_ods, t_ti_step_diag
   implicit none
   class(c_tint), intent(inout) :: ti   !< time integrator
   class(c_ode),  intent(in)    :: ode  !< ODE problem
   real(wp),      intent(in)    :: t    !< time level
   class(c_stv),  intent(in)    :: unow !< \f$u_n\f$
   class(c_ods),  intent(inout) :: ods  !< scratch (diagnostic vars.)
   class(c_stv),  intent(inout) :: unew !< \f$u_{n+1}\f$
   type(t_ti_step_diag), optional, intent(inout) :: step_diag
  end subroutine i_step
 end interface
 !> Generic clean-up subroutine of the time integrators
 abstract interface
  pure subroutine i_clean(ti , ode , step_diag)
   import :: c_tint, c_ode, t_ti_step_diag
   implicit none
   class(c_ode),  intent(in)    :: ode  !< ODE problem
   class(c_tint), intent(inout) :: ti   !< time integrator
   type(t_ti_step_diag), optional, intent(inout) :: step_diag
  end subroutine i_clean
 end interface

 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_time_integrators_base_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_time_integrators_base'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_time_integrators_base_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false. ) .or. &
          (mod_kinds_initialized.eqv..false. ) .or. &
      (mod_mpi_utils_initialized.eqv..false. ) .or. &
     (mod_state_vars_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_time_integrators_base_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_time_integrators_base_initialized = .true.
 end subroutine mod_time_integrators_base_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_time_integrators_base_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_time_integrators_base_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_time_integrators_base_initialized = .false.
 end subroutine mod_time_integrators_base_destructor

!-----------------------------------------------------------------------
 
 !> Override this subroutine if you use implicit integrators
 subroutine solve(x,ode,t,sigma,b,xl,ods)
  class(c_ode), intent(in)    :: ode
  real(wp),     intent(in)    :: t, sigma
  class(c_stv), intent(in)    :: b, xl
  class(c_ods), intent(inout) :: ods
  class(c_stv), intent(inout) :: x

  character(len=*), parameter :: &
    this_sub_name = 'solve'
  character(len=*), parameter :: &
    err_msg(4) = (/ &
   "Please, consider that you have to provide a real implementation  ",&
   "for this subroutine if you want to use implicit time integrators.",&
   "This subroutine in mod_time_integrators_base is provided only to ",&
   "simplify the implementation with explicit integrators.           " &
                 /)

   call error(this_sub_name,this_mod_name,err_msg)
 end subroutine solve
 
!-----------------------------------------------------------------------

end module mod_time_integrators_base

